A Molecular Dynamics Simulation of Heat Conduction of a Finite Length Single-walled Carbon Nanotube

Heat conduction of finite length single walled carbon nanotubes (SWNTs) was simulated by the molecular dynamics method with Tersoff-Brenner bond order potential. Temperature at each end of a SWNT was controlled by the phantom technique, and the thermal conductivity was calculated with Fourier’s law from the measured temperature gradient and the energy budgets in phantom molecules. The measured thermal conductivity did not converge to a finite value with increase in tube length up to 404 nm, but an interesting power law relation was observed. The phonon density of states and photon dispersion relations were directly extracted from simulation results for further analysis of heat conduction mechanism based on the phonon concept. INTRODUCTION The discovery of single-walled carbon nanotubes (SWNTs) in 1993 [1] has opened the fascinating new scientific and technological research fields. Many exciting physical properties of this new form of carbon have been revealed [2, 3]. Among them, the thermal conductivity of carbon nanotubes, which is speculated to be higher than any other material along the 1 Address correspondence to Shigeo Maruyama, Department of Mechanical Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan. Tel: +81-3-58416421, Fax: +81-3-5800-6983, E-Mail: [email protected] cylindrical axis [4], is examined in this paper. Recently, measurements of thermal conductivity of a 5 μm thick deposited “mat” of SWNTs were reported for randomly oriented [5] or magnetically aligned conditions [6]. Comparing with the temperature dependence of electrical conductance in the same condition, it was concluded that the contribution of electrons to the thermal conductivity is negligible in all temperature range [5, 6]. The direct measurements of individual nanotube (multi-walled carbon nanotubes, MWNTs) were also tried with the scanning thermal microscopy [7] and with MEMS assisted new measurement technology [8]. On the other hand, several molecular dynamics simulations [4, 9, 10] showed very high thermal conductivity such as 6600 W/mK at 300 K [4]. However, the estimated values of thermal conductivity were widely different from one another, even though the same Tersoff-Brenner potential [11] was employed. The most important difficulty was the convergence in frequency space and some extrapolations. Here, a finite length SWNT without the use of periodic boundary condition was simulated. The length-dependent thermal conductivity was calculated in the range of a few hundreds nanometers, which might be the size of future nanotube devices. In addition, the thermal conductivity of infinitely long nanotube could be approached with collections of physically meaningful length dependent results. The reported peaking phenomena [4, 9] of temperature dependence of thermal conductivity were also curious because no quantum corrections in the heat capacity were applied. Some simulation results about the temperature dependence are addressed. Another purpose of this study is the preliminary connection of molecular dynamics techniques to the solid-state heat conduction usually discussed as “phonon transport” in solid physics [12]. In principle, the molecular dynamics simulation can be used to obtain information for phonon transport dynamics such as phonon dispersion relation, group velocity, mean free path, boundary scattering rate and the rate of phonon-phonon scattering (Umklapp process). Especially, the phonon scatting at an interface is very important issue in recent thin film technology [13]. It is also anticipated that by developing the phonon concept to more general form in order to understand the thermal boundary resistance even in the liquid-solid interface [14]. As the first step, the direct calculation of phonon density of states and phonon dispersion relations from the molecular trajectories are demonstrated in this paper. SIMULATION TECHNIQUE Three armchair type SWNT structures with different diameters (5, 5), (8, 8), and (10, 10) in Figure 1 were chosen. Here, the chiral vector (n, m) (or chiral index) uniquely determines the geometrical structure of SWNTs [2, 3]. SWNTs with (10, 10) have been most widely studied because it was speculated that the laser-oven grown SWNTs [15] might have predominantly (10, 10) structure. However, recent studies showed the fairly random distribution of chiralities [16]. On the other hand, (5,5) has the almost same diameter as C60, and the large-scale production of SWNTs with this diameter is anticipated with the new generation technique using high-pressure and high-temperature CO gas [17]. The Brenner potential [11] with the simplified form [18] was employed as the potential function with the parameter set II (table 2 in [11]), which better reproduced the force constant. The velocity Verlet method was adopted to integrate the equation of motion with the time step of 0.5 fs. By applying the phantom heat bath model [19] to each end of a SWNT, the temperature difference was applied as in Figure 1(e). While one unit-cell molecules were fixed, one unit-cell molecules next to them were controlled by the Langevin equation. Here, no periodic boundary condition was applied to handle the finite size effect of carbon nanotubes. For the normal infinite length simulation, it is often discussed that the cell length of periodic boundary condition should be larger than the “mean free path” of phonon, which is argued to be order of 1 μm. Here, the physically meaningful nanotube length was varied from 6 nm through 404 nm for (5, 5) nanotube. With our configuration, thermal conductivity was calculated from the measured temperature gradient and the heat flux obtained by the energy budgets of phantom molecules. RESULTS AND DISCUSSIONS Thermal Conductivity An example of the temperature distribution along the tube is shown in Figure 2 for 202 nm long (5, 5) structure case (16000 atoms). Phantom temperatures at each end were set as 290 K and 310 K. After obtaining the average temperature of about 300 K with the auxiliary velocity scaling control, typically 1 ns simulations were performed for the equilibration with only phantom temperature control. Then, 1 ~ 2 ns calculation was used for the measurement of temperature distribution. The convergence of the temperature distribution was terribly poor with this technique especially for shorter nanotubes, so we disregarded the data for shorter nanotubes than 12.4 nm. With energy budgets of controlling phantom molecules, the heat flux along the tube can be simply calculated. Combined with the temperature gradient such as in Figure 2, the thermal conductivity λ can be calculated through Fourier’s equation: ) / ( z T q ∂ ∂ − = λ . As the cross-sectional area, a ring of van der Waals thickness 3.4 Å was used. The dependence of the thermal conductivity on the nanotube length is summarized in Figure 3. The thermal conductivity was diverging with the power-law characteristics with nanotube length [20] at least up to the 0.4 μm long nanotube for (5, 5). This divergence is most dominant for the smallest diameter nanotube (5, 5). This very striking behavior of thermal conductivity for (5, 5) is similar to the one-dimensional model calculations of thermal conductivity [21] where the divergence of λ with the power of 0.35 or 0.4 is discussed. It seems that the one-dimensional feature of heat conduction is really possible with the real material: the small diameter carbon nanotube. The thermal conductivity may converge when the tube length is much longer than the mean free path of energy carrying phonon. However, the thermal conductivity of about 600 W/mK for about 0.4 μm tube is still steadily increasing with the exponent of 0.27. The temperature dependence of measured thermal conductivity is plotted in Figure 4. In solid state physics textbook [12], thermal conductivity λ by the crystal vibration is expressed as ∑ = lv cv λ , where cv, l and v are heat capacity, mean free path, and group velocity of phonons, respectively. Usually, it is assumed that the temperature dependence of the group velocity or the sound velocity is negligible. At high temperature the phonon-phonon scattering known as Umklapp process is the dominant factor to reduce the phonon mean free path, and the thermal conductivity is expected to reduce as 1/T. Calculated results in Figure 4 approximately follow this relation, but even at very low temperature about 20 K without a peak. On the other hand, experimental results of MWNTs [8] or carbon fibers [22] showed the maximum thermal conductivity at around 300 K. At low temperature the phonon mean free path is expected to have almost constant value limited by the system size. On the other hand, the heat capacity varies as T at low temperatures because only long wavelength (low energy) phonons having T kB < ω h are excited to any appreciable extent at low temperature T. At high temperature, the heat capacity will approach to the constant value of 24.943 J/(mol K) approximately beyond the Debye temperature. Since the Debye temperature for carbon nanotube is expected to be as high as diamond, such as 2500 K, the quantum effect in heat capacity is very important even at room temperature. Since the classical molecular dynamics simulation cannot reproduce the quantum effect in heat capacity [23], the reported maximum value of thermal conductivity at 100 K [4] or around 300 K [9] are questionable. These coincidences are probably because the limit of the phonon mean free path due to the small system size. The relatively short nanotubes in Figure 4 also show the deviation from the 1/T relation. Phonon Density of States and Phonon Dispersion Relations As the first step to connect the molecular dynamics simulation to the phonon concept approach, phonon density of states and phonon dispersion relations were extracted from the simulated trajectories. The phonon density of states was measured as the power spectra of velocity fluctuations in Eq. (1). ( ) ) 0 ( ) ( exp ) ( α α α ∫ ω − = ω v t v t i dt D (1) Here, α takes r, θ, z for each velocity component in the cylindrical coordinate. The calculated density of states is shown in Figure 5 as the right hand side inserts of phonon dispersion relations. For this calculation, 8192 temporal data points saved every 5 fs were used. The photon dispersion relations were also directly measured as the time-space 2-D Fourier transforms in Eq. (2) of the displacement vector r’(z, t) from the equilibrium position re(z); ) ( ) , ( ) , ( ' z t z t z e r r r − = . ∫ ω − = ω α α ) exp( ) , ( ' ) , ( ' t i ikz t z r dt k R (2) Here, k is the wave vector along z-axis. Position vector r(z, t) of atoms in a line along zdirection separated by c c a a − = 3 (see Figure 1(d)) are used. In the case of Figure 5, 128 data points in z direction multiplied by 8192 temporal points were used. The phonon dispersion relation obtained from eigen values and eigen vectors of dynamical matrix made with the force-constant tensor scaled from 2 D graphite [24] is shown in Figure 5(d) as the comparison. The phonon modes from radial displacement or radial velocity component in Figure 5(a) have relatively low frequency including the radial breathing mode [25]. On the other hand, the vibrational modes in the graphitic plane appear as about 1600 cm peak in Figure 5(b) and (c). The unit cell of (5, 5) nanotube is made of 20 atoms (Figure 1(c)), and there are 60 vibrational modes. Group theory for this D5d symmetry structure shows that there are 12 non-degenerate and 24 doubly degenerate phonon branches [2]. Hence, there are distinguishable 36 branches of phonon dispersions in Figure 5(d). All of them were projected in Figure 5(a), (b) or (c). The group velocities of four acoustic phonons were estimated as longitudinal acoustic mode: 17 km/s (z-direction), transverse acoustic mode (degenerate): 7 km/s (r-direction), twisting acoustic mode: 10 km/s (θ-direction), in good agreement with the dynamical matrix result in Figure 5(d). SUMMARY Heat conduction of finite length single-walled carbon nanotubes (SWNTs) was simulated by the molecular dynamics method with the Tersoff-Brenner bond order potential. Without the use of periodic boundary conditions in the nanotube axis direction, thermal conductivity of finite length nanotube was calculated. The thermal conductivity was diverging with the power-law characteristics with nanotube length at least up to the 0.4 μm long nanotube for (5, 5). This feature can be compared with the theoretical 1 dimensional heat conduction calculations. The temperature dependency of measured thermal conductivity behaved as 1/T in high temperature range. The limitation of classical molecular dynamics in low temperature calculation was also discussed. As the first step for exploration the basic heat conduction mechanism, phonon density of states and phonon dispersion relations were directly extracted from the molecular dynamics trajectories. NOMENCLATURE a: Translation lattice vector ac-c: Bond length of carbon, nm c: velocity of light, m/s D(ω): Phonon density of states h : Planck constant, Js k: Wave vector, 1/nm kB: Boltzmann constant, J/K L: Length of a SWNT, nm q: Heat flux, W/m r: Cylindrical coordinate (radial direction), nm R(k,ω): Phonon dispersion relations T: Temperature, K t: Time, s v: Velocity, m/s x, y, z: Cartesian coordinates (z is along tube axis), nm Greek Symbols λ: Thermal conductivity, W/mK